#include "mhd.hh"
#include <map>
#include <math.h>

namespace Dune {
  namespace Equation {
    MHD_<1>::MHD_(Dune::ParameterTree& ptree){
      set_equ_paras(ptree);
    }
    MHD_<2>::MHD_(Dune::ParameterTree& ptree){
      set_equ_paras(ptree);
    }
    MHD_<3>::MHD_(Dune::ParameterTree& ptree){
      set_equ_paras(ptree);
    }
    MHD_<1>::~MHD_(){
    }
    MHD_<2>::~MHD_(){
    }
    MHD_<3>::~MHD_(){
    }

    void MHD_<1>::set_equ_paras(Dune::ParameterTree& ptree) {
      varname.push_back("density");
      varname.push_back("momentum1");
      varname.push_back("momentum2");
      varname.push_back("momentum3");
      varname.push_back("energy");  
      varname.push_back("b2");
      varname.push_back("b3");
      set_equ_para("gamma",5.0/3.0);
      set_equ_para("b1constant",0.0);  
    }
    void MHD_<2>::set_equ_paras(Dune::ParameterTree& ptree) {
      varname.push_back("density");
      varname.push_back("momentum1");
      varname.push_back("momentum2");
      varname.push_back("momentum3");
      varname.push_back("energy");  
      varname.push_back("b1");
      varname.push_back("b2");
      varname.push_back("b3");
      set_equ_para("gamma",5.0/3.0);
    }
    void MHD_<3>::set_equ_paras(Dune::ParameterTree& ptree) {
      varname.push_back("density");
      varname.push_back("momentum1");
      varname.push_back("momentum2");
      varname.push_back("momentum3");
      varname.push_back("energy");  
      varname.push_back("b1");
      varname.push_back("b2");
      varname.push_back("b3");
      set_equ_para("gamma",5.0/3.0);
    }

    void MHD_<1>::get_equation_flux(int whichdim, double* fequ, double* var){
      get_1d_fx(fequ,var);
    }
    void MHD_<2>::get_equation_flux(int whichdim, double* fequ, double* var){
      if(whichdim==0)
        get_2d_fx(fequ,var);
      if(whichdim==1)
        get_2d_fy(fequ,var);
    }
    void MHD_<3>::get_equation_flux(int whichdim, double* fequ, double* var){
      if(whichdim==0)
        get_2d_fx(fequ,var);
      if(whichdim==1)
        get_2d_fy(fequ,var);
      if(whichdim==2)
        get_3d_fz(fequ,var);  
    }

    double MHD<1>::get_thermal_pressure(double* var){
      double pressure,kinetic_e,magnetic_e;  
      kinetic_e=1.0/2.0*(std::pow(var[1],2)+std::pow(var[2],2)+std::pow(var[3],2))/var[0];
      magnetic_e=1.0/2.0*(std::pow(var[5],2)+std::pow(var[6],2)+std::pow(para["b1constant"],2));
      pressure=(para["gamma"]-1.0)*(var[4]-kinetic_e-magnetic_e);	
      return pressure;
    }
    double MHD<2>::get_thermal_pressure(double* var){
      double pressure,kinetic_e,magnetic_e;  
      kinetic_e=1.0/2.0*(std::pow(var[1],2)+std::pow(var[2],2)+std::pow(var[3],2))/var[0];
      magnetic_e=1.0/2.0*(std::pow(var[5],2)+std::pow(var[6],2)+std::pow(var[7],2));
      pressure=(para["gamma"]-1.0)*(var[4]-kinetic_e-magnetic_e);
      return pressure;
    }
    double MHD<3>::get_thermal_pressure(double* var){
      double pressure,kinetic_e,magnetic_e;  
      kinetic_e=1.0/2.0*(std::pow(var[1],2)+std::pow(var[2],2)+std::pow(var[3],2))/var[0];
      magnetic_e=1.0/2.0*(std::pow(var[5],2)+std::pow(var[6],2)+std::pow(var[7],2));
      pressure=(para["gamma"]-1.0)*(var[4]-kinetic_e-magnetic_e);
      return pressure;
    }

    double MHD<1>::local_max_eigenvalue(double* var){
      double tmp=0;
      double pressure=get_thermal_pressure(var);
      tmp = std::sqrt(para["gamma"]*pressure/var[0]
        +(pow(para["b1constant"],2)+pow(var[5],2)+pow(var[6],2))/var[0]);
      return tmp;
    }
    double MHD<2>::local_max_eigenvalue(double* var){
      double tmp=0;
      double pressure=get_thermal_pressure(var);
      tmp = std::sqrt(para["gamma"]*pressure/var[0]
        +(pow(var[5],2)+pow(var[6],2)+pow(var[7],2))/var[0]); 
      return tmp;
    }
    double MHD<3>::local_max_eigenvalue(double* var){
      double tmp=0;
      double pressure=get_thermal_pressure(var);
      tmp = std::sqrt(para["gamma"]*pressure/var[0]
        +(pow(var[5],2)+pow(var[6],2)+pow(var[7],2))/var[0]); 
      return tmp;
    }

    void MHDDecorate::get_1d_fx(double* fequ, double* var){
      double pressure,magnetic_e, bdotv;
      pressure = get_thermal_pressure(var);
      magnetic_e=1.0/2.0*(std::pow(var[5],2)+std::pow(var[6],2)+std::pow(para["b1constant"],2));  
      bdotv=(para["b1constant"]*var[1]+var[5]*var[2]+var[6]*var[3])/var[0];
      fequ[0]=var[1];
      fequ[1]=std::pow(var[1],2)/var[0]+pressure+magnetic_e-std::pow(para["b1constant"],2.0);
      fequ[2]=var[1]*var[2]/var[0]-para["b1constant"]*var[5];
      fequ[3]=var[1]*var[3]/var[0]-para["b1constant"]*var[6];
      fequ[4]=(var[4]+pressure+magnetic_e)*var[1]/var[0]-para["b1constant"]*bdotv;
      fequ[5]=(var[5]*var[1]-para["b1constant"]*var[2])/var[0];
      fequ[6]=(var[6]*var[1]-para["b1constant"]*var[3])/var[0];  
    }

    void MHDDecorate::get_2d_fx(double* fequ, double* var){
      double pressure,magnetic_e, bdotv;
      pressure = get_thermal_pressure(var);
      magnetic_e=1.0/2.0*(std::pow(var[5],2)+std::pow(var[6],2)+std::pow(var[7],2));
      bdotv=(var[5]*var[1]+var[6]*var[2]+var[7]*var[3])/var[0];
      fequ[0]=var[1];
      fequ[1]=std::pow(var[1],2)/var[0]+pressure+magnetic_e-std::pow(var[5],2.0);
      fequ[2]=var[1]*var[2]/var[0]-var[5]*var[6];
      fequ[3]=var[1]*var[3]/var[0]-var[5]*var[7];    
      fequ[4]=(var[4]+pressure+magnetic_e)*var[1]/var[0]-var[5]*bdotv;
      fequ[5]=0.0;
      fequ[6]=(var[6]*var[1]-var[5]*var[2])/var[0];
      fequ[7]=(var[7]*var[1]-var[5]*var[3])/var[0];
    }

    void MHDDecorate::get_2d_fy(double* fequ, double* var){
      double pressure,magnetic_e, bdotv;
      pressure = get_thermal_pressure(var);
      magnetic_e=1.0/2.0*(std::pow(var[5],2)+std::pow(var[6],2)+std::pow(var[7],2));
      bdotv=(var[5]*var[1]+var[6]*var[2]+var[7]*var[3])/var[0];    
      fequ[0]=var[2];
      fequ[1]=var[2]*var[1]/var[0]-var[6]*var[5];
      fequ[2]=std::pow(var[2],2)/var[0]+pressure+magnetic_e-std::pow(var[6],2.0);        
      fequ[3]=var[2]*var[3]/var[0]-var[6]*var[7];
      fequ[4]=(var[4]+pressure+magnetic_e)*var[2]/var[0]-var[6]*bdotv;
      fequ[5]=(var[5]*var[2]-var[6]*var[1])/var[0];
      fequ[6]=0.0;    
      fequ[7]=(var[7]*var[2]-var[6]*var[3])/var[0];
    }

    void MHDDecorate::get_3d_fz(double* fequ, double* var){
      double pressure,magnetic_e, bdotv;
      pressure = get_thermal_pressure(var);
      magnetic_e=1.0/2.0*(std::pow(var[5],2)+std::pow(var[6],2)+std::pow(var[7],2));
      bdotv=(var[5]*var[1]+var[6]*var[2]+var[7]*var[3])/var[0];    
      fequ[0]=var[3];
      fequ[1]=var[3]*var[1]/var[0]-var[7]*var[5];
      fequ[2]=var[3]*var[2]/var[0]-var[7]*var[6];    
      fequ[3]=std::pow(var[3],2)/var[0]+pressure+magnetic_e-std::pow(var[7],2.0);            
      fequ[4]=(var[4]+pressure+magnetic_e)*var[3]/var[0]-var[7]*bdotv;
      fequ[5]=(var[5]*var[3]-var[7]*var[1])/var[0];
      fequ[6]=(var[6]*var[3]-var[7]*var[2])/var[0];
      fequ[7]=0.0;    
    }
  } // Equation
} // Dune
